clc;clear all;close all;
load('solN2_ded.mat')

%ALPHall=[1 0.8 0.6 0.5 0.45 0.4 0.35 0.3 0.29 0.28 0.27];

ALPHall2=[1 0.7 0.5 0.3 0.2 0.1 0.05 0.03 0.02 0.012 0.005 0.002 0]

ALPHA=1;
X=XX_ded(1,:);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
for i=1:T;
   X(T+i) = L2(lamda,X(T+i),phi(2)); %l    
end
X(3*T+1)=-X(3*T+1);
X(3*T+3+N-2:3*T+3+2*(N-2))=1/lamda;

kistart = flip(kj);
phi=flip(phi);

for ii=1:length(ALPHall)

ALPHA=ALPHall2(ii)
if (ii>1)
    X=XX_ded2(ii-1,:);
end
X(3*T+3+N-2:3*T+3+2*(N-2))
X(3*T+1)
if ii==2
    X(3*T+3+N-2:3*T+3+2*(N-2))=1.6; %lambda
end
if ii==3
    X(3*T+3+N-2:3*T+3+2*(N-2))=1.65; %lambda
end
if ii==4
    X(3*T+3+N-2:3*T+3+2*(N-2))=1.5; %lambda
    X(3*T+1)=-0.5; %Delta1
end
if ii==5
    X(3*T+3+N-2:3*T+3+2*(N-2))=1.6; %lambda
    X(3*T+1)=-0.4; %Delta1
end
if ii==6
    X(3*T+1)=-0.57; %Delta1
end
if ii==7 | ii==8
    X(3*T+3+N-2:3*T+3+2*(N-2))=1.56; %lambda
    X(3*T+1)=-0.4; %Delta1
end
if ii==9 | ii==11 | ii==12
    X(3*T+3+N-2:3*T+3+2*(N-2))=1.57; %lambda
    X(3*T+1)=-0.4; %Delta1
end
if ii==13
    X(3*T+3+N-2:3*T+3+2*(N-2))=1.57; %lambda
    X(3*T+1)=-0.3; %Delta1
end

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
X1=X;
sum(resid3v(X1).^2)

X=X1;
tolf0 = 1e-7;
[X,check] = broydn('resid3v',X,tolf0)%,0,1);
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
X2=X;

sumresid_all2(ii)=sum(resid3v(X).^2)

XX_ded2(ii,:)=X;
    
k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = -Delta1
deductible = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
deductible_all2(ii)=deductible;

kstst = fzero('findstst',kstart+0.3,options);

rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end
wstst = (1-alpha)*kstst^alpha*estst^(-alpha); %w=F_e
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst.^alpha*estst.^(1-alpha) - depr*kstst - govexp); %c from resource constraint

klast = [kstart X(1:T-1)]; %k_{t-1}
e=[l' (l'*(lamda.^(-sigmaC/sigmaL).*(phi(2:N)/phi(1))'.^(1/sigmaL)))]*(popweight'.*phi);
e=e';

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];
lnext = [X(T+2:2*T) lstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
listst = [lstst, L2(lamda,lstst,phi(2:N))']; %hours worked of groups in the long run
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt), L2(lamda,l(tt),phi(2:N))']; %hours worked of groups during the transition
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)- omega*li.^(1+sigmaL)/(1+sigmaL)))...
    + betta^(T)/(1-betta)*(log(cistst) - omega*listst.^(1+sigmaL)/(1+sigmaL));
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC) - omega*li.^(1+sigmaL)/(1+sigmaL)))...
    + betta^(T)/(1-betta)*((cistst.^(1-sigmaC)-1)/(1-sigmaC) - omega*listst.^(1+sigmaL)/(1+sigmaL));
end

if sum(V>=Vsq)==2
    PI(ii)=1
else
    PI(ii)=0
end
Vall_ded2(ii,:)=V
Vsq

%tax rates
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
plot(tauKt)
[swi dur]=min((tauKt-taumax/2).^2);
duration2(ii)=dur
deductible_all2

save('solN2PO_ded.mat')

end

Vcap=Vall_ded(:,1);
Vwor=Vall_ded(:,2);

term1=(1-sigmaC)*((1-betta)*Vcap+omega*listst_sq(1)^(1+sigmaL)/(1+sigmaL));
epscap=((term1+1).^(1/(1-sigmaC))/cistst_sq(1)-1)*100;
term2=(1-sigmaC)*((1-betta)*Vwor+omega*listst_sq(2)^(1+sigmaL)/(1+sigmaL));
epswor=((term2+1).^(1/(1-sigmaC))/cistst_sq(2)-1)*100;

Vcap=Vall_ded2(:,2);
Vwor=Vall_ded2(:,1);

term1=(1-sigmaC)*((1-betta)*Vcap+omega*listst_sq(1)^(1+sigmaL)/(1+sigmaL));
epscap2=((term1+1).^(1/(1-sigmaC))/cistst_sq(1)-1)*100;
term2=(1-sigmaC)*((1-betta)*Vwor+omega*listst_sq(2)^(1+sigmaL)/(1+sigmaL));
epswor2=((term2+1).^(1/(1-sigmaC))/cistst_sq(2)-1)*100;

epscap_ded=[flip(epscap2)' epscap']
epswor_ded=[flip(epswor2)' epswor']

save('solN2PO_ded.mat')

clc;clear all;close all;
load('solN2PO_ded_taumax.mat')
load('solN2PO_ded.mat')

plot(epswor_ded,epscap_ded,'-k','LineWidth',2.5);
hold on
plot(epswor_ded_taumax,epscap_ded_taumax,'-.k','LineWidth',2.5);
xline(0,':k','LineWidth',1.5)
yline(0,':k','LineWidth',1.5)
xlabel("workers' welfare increase (percent)",'FontSize',12)
ylabel("capitalists' welfare increase (percent)",'FontSize',12)
legend('PO with (negative) deductible','PO with (negative) deductible and $\tau^k_t=\tilde{\tau}, \forall t$','Interpreter','latex')
